Modelling of transient wind turbine loads during pitch motion

In connection with the design of wind turbines and their control algorithms, the transient loads, especially generated by time varying blade loads, are very important. In the Blade Element Momentum method, the most widespread tool in the wind turbine industry, the time constants necessary to describe these problems are not an inherent part of the model. In the present study two different approaches are used to determine these time constants, namely the EllipSys3D Navier-Stokes solver, and an actuator disc method. The time constants estimated by the actuator disc method is afterwards used in a standard BEM method and a BEM method coupled with a near wake method. The resulting transients are compared to measured values, for a step pitch change taken from the unique data set of the NREL/NASA Ames experiment, with good agreement for the case with low loading. Following this, a series of pitch steps corresponding to more realistic operational conditions are investigated. Introduction Dealing with transient aerodynamic loads for wind turbine these are typically handled at two levels, unsteady airfoil aerodynamics, and the so called dynamic inflow models. The first part covers the unsteady non-separated effects from shed vorticity [1,2] and the dynamic stall models, phenomena which both deals with relatively fast motion with a time scale in the order of the chord divided by the velocity seen by the airfoil. The second part of the dynamic loads, the dynamic inflow, or dynamic induction, is phenomena related to variations in the inflow velocity, structural vibrations and tip pitch changes. The time scales important for this are of size rotor radius divided by the free stream velocity. In the present paper the main focus is on the second part of these models namely the dynamic inflow models. Dynamic inflow models have been developed early on in connection with helicopter aerodynamics for forward flight, see [3],[4]. It is important though to be aware that as stated in earlier work dealing with dynamic inflow, [5], that the wind turbine operational conditions are different from the conditions experienced in connection with helicopters. First of all, wind turbines are designed for power extraction at high induction factors ~1/3, while helicopters are designed for maximum thrust with minimum power consumption at low induction factors. As the dynamic inflow problem is closely related to the induction, and is most important for high induction factors, some of the assumptions used in connection with the helicopter related models may not be applicable for wind turbine aerodynamics. Basically, the dynamic inflow models describe the effect that when the load changes, for example due to a change in blade pitch, this will not be instantaneously reflected in the induced velocity field, but gradually a new equilibrium will be established between the load on the blade, the rotor wake and the induced velocity at the rotor plane. The direct consequence of the load overshoots experienced in connection with fast pitch changes, are of minor importance for the necessary structural strength of the turbine and the extreme loads. The main area of concern for these phenomena is the influence on the damping properties in connection with dynamic instabilities and vibrations. The importance of the values of the damping coefficients was investigated in [6], showing that neglecting the transient behavior will result in overestimating of the aerodynamic damping, that may result in unwanted vibrations in the real turbine. When investigating problems related to dynamic inflow, one problem is to obtain accurate and detailed measurements. One source for this has been the Tjæreborg wind turbine experiment [7, 8, 9, 10, 11]. Even though these measurements have very good quality, they lack detailed information about the spanwise load distribution along the blade. In connection with the NREL/NASA Ames wind tunnel test [12, 13], detailed measurements were taken for a series of operational conditions [14] ranging from upwind axial operation, down wind operation, yaw operation and a series of step pitch operation that are interesting in connection with the present investigation. In the current work, one series of step pitch changes taken from [15] are investigated using both a full 3D CFD code, and two versions of the BEM model, with the main focus of looking at the time constants for the settling of the loads. Following the comparison with measurements, a parametric study is performed to investigate the dynamic inflow effects for other operational conditions. Code description/Method The in-house flow solver EllipSys3D used in the present study for the CFD computations, is developed in cooperation between the Department of Mechanical Engineering at DTU and The Department of Wind Energy at Risø National Laboratory, see [16, 17 and 18]. The EllipSys3D code is a multi-block finite volume discretization of the incompressible Reynolds Averaged Navier-Stokes (RANS) equations in general curvilinear coordinates, and is second order accurate in both time and space. The code is parallelized with MPI for executions on distributed memory machines, using a non-overlapping domain decomposition technique. In the present work the turbulence in the boundary layer is modeled by the k-ω SST eddy viscosity model [19]. The rotation and pitch of the rotor is modeled using a moving mesh formulation in a fixed frame of reference. The moving mesh option has been implemented in the EllipSys3D solver in a generalized way allowing arbitrary deformation of the computational mesh, following [20]. It has been used for doing unsteady simulations for several years both for stiff rotors in yaw and fully coupled aeroelastic computations [21, 22, 23, 24]. The mesh used in the present study has 5.2 million points and has previously been used in connection with yaw computations [21], where more details can be found. Additionally, two engineering models are used in the present study; a Blade Element Momentum model (BEM model) and a BEM model coupled to a near-wake model (NW model). The BEM model is well-known and is used in almost all aeroelastic models to compute the induction, [25]. Basically, the BEM model is a steady model and therefore a sub-model must be introduced to compute the unsteady induction caused by load changes on the rotor due to e.g. pitch changes or eigen-motion of the rotor or the blades. This submodel must simulate the delay in changes of induction at the rotor disc when going from one loading on the disc to another loading. Such a load change will cause a change in trailed and shed vorticity from the rotor blades and the new state is first stationary when the whole wake system of the rotor corresponds completely to the new loading on the rotor. In the present implementation of the dynamic induction model a filter function is applied on the instantaneous computed induced axial velocity at the rotor disc. Both a 1 order and a 2 filter have been tested. The filter is modeled using the indicial function technique. Besides the model for filtering the induced velocity it is necessary with a sub-model for computating the effect of the shed vorticity in order to model the complete influence of the vortex system in the wake properly. This is also done with an indicial function technique and follows the implementation in the Beddoes-Leishman dynamic stall model, [26]. The near wake model was originally developed by Beddoes, [27], for computation of high time resolution air loads on helicopter rotors, and is a way to include the radial dependency of the loads. Recently it has been implemented at Risø for use on wind turbine rotors, [28]. In short, the main idea with the model is to simulate the downwash from the first 90 degree of the trailed vortex system behind the individual blades with a simple, unsteady lifting line model. Again the indicial function technique is used to compute the downwash in the model. As the near wake model only computes the induction from the first part of the trailed vortex system, a model to compute the induction from the far wake is needed. In the present case a BEM model is used but where the loading in the form of the local thrust coefficient is scaled down with a constant which in the present case is 0.85. This constant was found correlating the load distribution in the combined NW model with the result of a standard BEM computation. It should be noted that no tip correction is used in the NW model as the trailed vortex system of the individual blades are modeled and forces the loading at the tip to approach zero. The same model for the shed vorticity as described above is also used in the NW model. Actuator disc simulations with a uniform loading and with a step change from one thrust coefficient to another has been used as a basis for estimation of time constants to be used in both engineering models. For the NW model only one time constant is used as the NW model itself takes care of the initial fast decay, while the slower decay must be handled by the filtering of the induced velocities. Normally, the time constants are made dimensionless using the free stream velocity and the radius of the rotor. Another option used in the present work is to use the instantaneous velocity in the wake for the normalization R a V ) 5 . 1 1 ( ~ − = ∞ τ τ . In both the BEM and the NW model, this is the approach used, and as will be seen later in the paper, this assures that the time constants reflects the slower development of the wake when high induction are present. The effect of changing the normalization can be seen in Figure 1, where the use of the actual wake velocity makes the curves collapse. Figure 1: The dependency of the non-dimensionale time constant on the choice of normalization velocity used, left figure using the free stream velocity, right figure shows the use of the actual wake velocity. In the present work, when analyzing the transient behavior of the dynamic inflow effects of both the experiment and the numerical simulations, the dynamic inflow effects are parameterized as a sum of two exponential decay functions, as shown below. The sum of two exponential functions is used, instead of just a single one, as a compromise to allow us to capture both the fast decay initially and the more slow decay in the last stage of the development towards equilibrium. ) / ) ( exp( ) / ) ( exp( ) 1 ( 2 0 1 0 τ α τ α χ t t t t − − + − − − = Alternatively an expression where the time ‘constant’ is directly a function of time could be used. The normal forces, Fn, are normalized using the following expression: